Single-cell profiling of response to neoadjuvant chemo-immunotherapy in surgically resectable esophageal squamous cell carcinoma

Background The efficacy of neoadjuvant chemo-immunotherapy (NAT) in esophageal squamous cell carcinoma (ESCC) is challenged by the intricate interplay within the tumor microenvironment (TME). Unveiling the immune landscape of ESCC in the context of NAT could shed light on heterogeneity and optimize therapeutic strategies for patients. Methods We analyzed single cells from 22 baseline and 24 post-NAT treatment samples of stage II/III ESCC patients to explore the association between the immune landscape and pathological response to neoadjuvant anti-PD-1 combination therapy, including pathological complete response (pCR), major pathological response (MPR), and incomplete pathological response (IPR). Results Single-cell profiling identified 14 major cell subsets of cancer, immune, and stromal cells. Trajectory analysis unveiled an interesting link between cancer cell differentiation and pathological response to NAT. ESCC tumors enriched with less differentiated cancer cells exhibited a potentially favorable pathological response to NAT, while tumors enriched with clusters of more differentiated cancer cells may resist treatment. Deconvolution of transcriptomes in pre-treatment tumors identified gene signatures in response to NAT contributed by specific immune cell populations. Upregulated genes associated with better pathological responses in CD8 + effector T cells primarily involved interferon-gamma (IFNγ) signaling, neutrophil degranulation, and negative regulation of the T cell apoptotic process, whereas downregulated genes were dominated by those in the immune response-activating cell surface receptor signaling pathway. Natural killer cells in pre-treatment tumors from pCR patients showed a similar upregulation of gene expression in response to IFNγ but a downregulation of genes in the neutrophil-mediated immunity pathways. A decreased cellular contexture of regulatory T cells in ESCC TME indicated a potentially favorable pathological response to NAT. Cell–cell communication analysis revealed extensive interactions between CCL5 and its receptor CCR5 in various immune cells of baseline pCR tumors. Immune checkpoint interaction pairs, including CTLA4-CD86, TIGIT-PVR, LGALS9-HAVCR2, and TNFSF4-TNFRSF4, might serve as additional therapeutic targets for ICI therapy in ESCC. Conclusions This pioneering study unveiled an intriguing association between cancer cell differentiation and pathological response in esophageal cancer patients, revealing distinct subgroups of tumors for which neoadjuvant chemo-immunotherapy might be effective. We also delineated the immune landscape of ESCC tumors in the context of clinical response to NAT, which provides clinical insights for better understanding how patients respond to the treatment and further identifying novel therapeutic targets for ESCC patients in the future. Supplementary Information The online version contains supplementary material available at 10.1186/s13073-024-01320-9.


Background
Esophageal cancer (EC) is the eighth most frequent cancer and the sixth leading cause of cancer-related mortality worldwide [1].EC has two predominant histological subtypes, esophageal adenocarcinoma (EAC), and esophageal squamous cell carcinoma (ESCC), with distinct geographical and racial variabilities.Unlike Western nations, ESCC is particularly more prevalent in Central and Eastern China, constituting approximately 90% of all EC cases [2].Despite incremental advances in diagnostics and therapeutics, ESCC has a dismal 5-year survival rate of 12-20% [3].Neoadjuvant chemoradiotherapy followed by surgery has become the standard treatment for patients with resectable locally advanced ESCC [4].However, nearly half of patients develop local recurrence or distant metastases following surgery, with an anticipated increase in toxicity levels and severe side effects [5,6].Therefore, there is an unmet clinical need to explore novel and effective treatments to improve patient survival.
In recent years, incorporating immune checkpoint inhibitors (ICIs) that target programmed cell death 1 (PD-1) with its ligand PD-L1 has emerged as a promising neoadjuvant treatment strategy in early-stage solid tumors, such as breast cancer, lung cancer, and ESCC [7][8][9].Several clinical trials, including KEYNOTE 590 [10] and CheckMate 649 [11], have demonstrated promising antitumor activity and the safety of immunotherapies with or without chemotherapy in patients with advanced ESCC.Nevertheless, despite the encouraging preliminary results yielded by the use of ICIs in preoperative neoadjuvant therapy settings, several unsolved issues require attention.The identification of biomarkers that reflect intricate interactions between the tumor and the immune system will be pivotal in effectively stratifying patients who will truly benefit from ICI neoadjuvant therapy.
Positive PD-L1 expression has been reported in 18.9 to 45% of ESCC patients [12][13][14][15].Although the relationship between PD-L1 expression and clinical outcomes in ESCC has been previously assessed, there is ongoing debate regarding whether PD-L1 expression is positively or negatively correlated with prognosis [16][17][18].Theoretically, the extent of CD8 + T cell infiltration within the TME is associated with ICI response, given that blocking the PD-1/PD-L1 axis can rejuvenate the cytotoxic effect of CD8 + exhausted T cells [19].However, recent studies harnessing neoadjuvant immunotherapy showed no significant difference in CD8 + T cells between responders and nonresponders [20][21][22][23].Furthermore, prior studies have demonstrated that the upregulation of inflammatory and interferon-gamma (IFNγ) signaling-related gene signatures in inoperable EAC patients serve as on-treatment markers of ICI efficacy, while the elevated expression of E2F targets and genes related to extracellular matrix has been associated with tumor growth and resistance to ICI treatment [24].The complex nature of IFNγ in the TME is continuously unfolding, and increasing evidence suggests that IFNγ may exert different functions depending on the immune context in which it is produced [25][26][27][28][29].
Single-cell RNA sequencing (scRNA-seq) has been arising as an unbiased, cutting-edge method that allows the grouping of cells based on their distinct transcriptional signatures without prior knowledge of genes and proteins of interest.Single-cell transcriptomic analysis offers an extraordinary opportunity to conduct an indepth analysis of heterogeneous cells, including malignant cells, immune cells, and stromal cells.Indeed, much attention is currently focused on deconvolving intertumoral and intra-tumoral heterogeneity, which might imply various mechanisms to evade antitumor immune responses.Several recent studies have extensively investigated the intricate TME in esophageal cancer patients at the single-cell level in the context of neoadjuvant therapies; however, these studies have been constrained by relatively small sample sizes [24,30,31].In this study, we undertake the scRNA-seq assay to uncover the key role immunity pathways.A decreased cellular contexture of regulatory T cells in ESCC TME indicated a potentially favorable pathological response to NAT.Cell-cell communication analysis revealed extensive interactions between CCL5 and its receptor CCR5 in various immune cells of baseline pCR tumors.Immune checkpoint interaction pairs, including CTLA4-CD86, TIGIT-PVR, LGALS9-HAVCR2, and TNFSF4-TNFRSF4, might serve as additional therapeutic targets for ICI therapy in ESCC.
of intra-tumoral cell type composition and their relationship to differential pathological responses to neoadjuvant therapy with PD-1 monoclonal antibodies plus chemotherapy in ESCC.

Participants and sample collection
Participants who had received immune therapy, such as monoclonal antibodies against programmed death-1 (PD-1)/PD-1 ligand (PD-L1) and cytotoxic T lymphocyte-associated antigen-4 (CTLA-4) or targeted therapies and radiotherapy were ineligible for this study.After screening, 22 patients who were histologically diagnosed with stage II to IV esophageal cancer between July 2020 and March 2021 and had not received the treatment mentioned above were included in this study.Relevant clinical information, such as sex, age, and treatment information, was retrieved from the database for each patient.The clinical stage was determined according to the American Joint Committee on Cancer (eighth edition) [32].This study was conducted following the Declaration of Helsinki and was approved by the Institutional Ethics Committee of Tangdu Hospital, Fourth Military Medical University (No. 202102-23).Written informed consent was obtained from each patient before sample collection.
All 22 patients received neoadjuvant chemo-immunotherapy (NAT) comprising tislelizumab (BeiGene Co. Ltd, Shanghai, China) (N = 20) or camrelizumab (Jiangsu Hengrui Pharmaceuticals Co. Ltd) (N = 2), along with carboplatin/nedaplatin (Qilu Pharmaceutical Co. Ltd, Jinan, China) and albumin-bound paclitaxel (Qilu Pharmaceutical Co. Ltd, Jinan, China) for three cycles (3 weeks/cycle).Among these, 18 patients underwent curative-intent complete surgery following NAT.Tissue samples were obtained by endoscopic biopsy under general anesthetic during routine staging (pre-NAT T_B and N_B) or sampling of esophagectomy resection (post-NAT T_A and N_A) [30,31].Adjacent normal esophageal tissue biopsies (N_B) or tissue specimens (N_A) were obtained at least 2.0 cm distant from the matched tumor (Fig. 1A).Pathological response to NAT was evaluated based on the presence of viable residual tumor cells (VRTCs) in the resected primary tumor [33].In particular, the presence of VRTCs ≥ 1% and ≤ 10% in the resected tumor and all resected lymph nodes was defined as the major pathological response (MPR), whereas the presence of ≤ 1% VRTCs in the resected specimen was defined as pathological complete response (pCR).The presence of > 10% VRTCs in the resected specimen was defined as an incomplete pathological response (IPR).Consequently, the study cohort was categorized into three groups: 7 cases of pCR, 6 cases of MPR, and 5 cases of IPR.This comprised a total of 32 tumor samples (16 T_B and 16 T_A, respectively) and 24 adjacent samples (11 N_B and 13 N_A, respectively) (Fig. 1A).After excluding poor-quality samples, there were 15 T_B samples, 7 N_B samples, 12 T_A samples, and 12 N_A samples subjected to scRNA-seq analysis (Additional file 1: Fig. S1A).Available samples were aggregated for uniform manifold projection (UMAP) dimensionality reduction to attain a holistic overview of the single-cell atlas of the ESCC TME.Principal findings regarding cancer and specific immune cell populations were exclusively based on pre-treatment tumor samples (T_B).

Library preparation and single-cell RNA sequencing
Freshly excised samples were stored in the sCelLiVE ® Tissue Preservation Solution (Singleron Biotechnology, Nanjing, China) on ice and transferred to the laboratory within 30 min after resection.All specimens were washed with Hanks Balanced Salt Solution (HBSS) three times, minced into 1-mm cubic pieces, and then digested with 3 mL sCelLiVE ® Tissue Dissociation Solution (Singleron Biotechnology, Nanjing, China) at 37 °C for 15 min.The cell suspension was collected and filtered through a 40-μm sterile strainer.To remove red blood cells, the GEXSCOPE ® red blood cell lysis buffer (RCLB) (Singleron Biotechnology, Nanjing, China) was added into the cell suspension at a ratio of 2:1 (RCLB: cell) and incubated at room temperature for 5 min.The mixture was then centrifuged at 300 × g for 5 min at 4℃.The supernatant was removed and the sediment was gently resuspended in phosphate-buffered saline (PBS) (HyClone, China).
scRNA-seq libraries were constructed as previously described [34].In brief, single-cell suspensions (~ 2 × 10 5 cells/mL) were loaded onto a microfluidic chip of the Singleron Matrix ® Single Cell Processing System (Singleron Biotechnology, Nanjing, China).After the mRNA from each cell was labeled with a unique molecular identifier (UMI), single-cell suspensions were collected from the microfluidic chip, followed by reverse transcription and PCR amplification to obtain the scRNA-seq libraries.Individual libraries were diluted to 4 nM and paired-end sequenced on the NovaSeq 6000 sequencing system (Illumina, San Diego, USA).

scRNA-seq data processing
The CeleScope (version 1.9.0, https:// github.com/ singl eron-RD/ CeleS cope) pipeline was used to generate the gene expression matrix containing normalized gene counts versus cells per sample.The expression matrix was filtered for each sample dataset through the Seurat (v3.1.2)R toolkit [35].The primary exclusion criteria are as follows: (1) cells with a gene count less than 200 or with a top 2% gene count; (2) cells with a maximum of 2% of UMI counts; (3) cells with mitochondrial content more than 50%; (4) genes expressed in less than five cells.The quality control data for scRNA-seq can be found in Additional file 2.

Dimensionality reduction, clustering, and cell type annotation
The filtered raw count matrix was normalized by dividing by the total counts per cell and transformed by calculating the natural logarithm.The top 2000 variable genes were selected using the "FindVariableFeautres" function in Seurat (v3.1.2) [35].Principal component analysis was then performed, and the top 20 principal components were subjected to dimensionality reduction and subsequent cell type clustering using the "FindClusters" function.Batch effect correction was performed using Harmony v1.0 [36].Cell clustering was achieved by the Louvain algorithm, setting the resolution parameter to 0.8, and visualized by the Uniform Manifold Approximation and Projection (UMAP) method [37].Cell-ID (v0.1.0,https:// github.com/ Rause llLab/ CelliD) [38] was used for automated annotation of cell types, followed by manual correction of cell types by classical markers in the SynEcoSys database (Singleron Biotechnology, Nanjing, China) [39][40][41].All Seurat functions were run with default parameters unless otherwise specified.

Inference of copy number variations from scRNA-seq data
Inferred copy number variations (CNVs) were analyzed using the inferCNV tool (v1.2.1; https:// github.com/ broad insti tute/ infer CNV) [42] using non-malignant cells as baselines.In brief, genes expressed in more than three cells were sorted by their genomic locations on each chromosome.A sliding window comprising 101 genes was used to smooth the relative expression on each chromosome to remove the effect of gene-specific expression.The relative expression values were centered at 1, and the ceiling of the relative expression values was set to 1.5 standard deviations from the residual normalized expression levels.The R package Pheatmap (v1.0.12; https:// github.com/ raivo kolde/ pheat map) [43] was used to visualize patients' inferred CNV profiles.

Hotspot analysis
Hotspot (v0.9.1, https:// github.com/ Yosef Lab/ Hotsp ot) was used to identify functional gene modules that illustrate heterogeneity within cancer cell subsets [46,47].Briefly, we used the "danb" model and selected the top 500 genes with the highest autocorrelation z-score for module identification.Modules were then identified using the "create_modules" function with the parameters: min_gene_threshold = 15 and fdr_threshold = 0.05.Module scores were calculated by using the "calculate_mod-ule_scores" function.Jaccard similarity coefficients were used to evaluate the transcriptional similarity between 12 cancer modules and signatures of cancer cell clusters.

SCENIC analysis
Single-cell regulatory network inference and clustering (SCENIC), a computational method for network inference and motif discovery allowing high-confidence prediction of key regulators and their direct target genes, was performed using the R package "SCENIC" (v0.1.5;https:// github.com/ aerts lab/ SCENIC) [48,49].The coexpression modules were run by GRNBoost.We downloaded the motifs database for Homo sapiens from https:// pysce nic.readt hedocs.io/ en/ latest/.The input matrix consisted of the normalized expression values of the cells of interest.

Differentially expressed gene analysis
Differentially expressed genes (DEGs) were identified by the "FindMarker" function in Seurat based on the Wilcoxon rank sum test.Genes expressed in more than 10% of cells in a cluster and with an average log (fold change) greater than 0.25 and an adjusted P value less than 0.05 were identified as DEGs.The cutoff for DEGs remains the same unless otherwise specified.

Pathway enrichment analysis
Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis were performed using the "clusterProfiler" package in R [50].GO gene sets include three categories, namely biological process (BP), cellular component (CC), and molecular function (MF).Pathways with an adjusted P value less than 0.05 were considered significantly enriched.

Cell communication analysis
The CellPhoneDB (v2.1.0),equipped with a known ligand-receptor pair database, was used to predict cellcell interaction [51,52].The number of ineffective distributions of reciprocal exchanges was set to 1000 and otherwise followed the default settings of the software.The predicted interaction pairs with a P value less than 0.05 and average log (fold change) greater than 0.1 were considered significant.

Survival analysis
The Cancer Genome Atlas (TCGA)-esophageal carcinoma (ESCA) dataset, which includes bulk RNA-seq data and overall survival information for 159 patients, was obtained from cBioPortal (https:// www.cbiop ortal.org/).To evaluate the relationship between the target gene set and patient prognosis, single sample gene set enrichment (ssGSVA) was used to calculate the score of each gene set for samples as previously described [53].Samples were assigned to either high or low-expression groups using the best score as the cutoff.Kaplan-Meier curves were used to show differences in the overall survival of patients between high and low-expression groups for each gene set using the "Survival" package in R (https:// github.com/ thern eau/ survi val).

Statistical analysis
The unpaired two-sample Wilcoxon test was used to compare cell distribution within two groups of cells.An unpaired two-sample t-test was used to compare the mean gene expression or gene signatures within two groups of cells.The Jonckheere-Terpstra test was used to assess the statistical significance of the trends across IPR, MPR, and pCR groups.Log-rank test was used to compare the overall survival of subgroup patients.A two-sided P value of less than 0.05 was considered significant for all tests unless indicated otherwise (*P < 0.05, **P < 0.01, ***P < 0.001).All statistical analyses were performed in R (version 4.1.3).

Patient overview
This study enrolled 22 stage II-IV ESCC patients, of whom 20 underwent neoadjuvant tislelizumab combined with chemotherapy, and two were treated with camrelizumab plus chemotherapy (Fig. 1A).The baseline clinical characteristics of patients are summarized in Table 1.In particular, 14 patients were male, and half were above 65 years of age.Only one patient was diagnosed with esophageal adenosquamous carcinoma (EASC).Among the 22 patients, 15 (68.2%) had PD-L1 TPS scores < 1%, 1 (4.5%) had TPS ranging from 1 to 50%, 4 (18.2%) had TPS scores ≥ 50%, and TPS scores of two patients were unknown.
Notably, four patients had an unknown pathological response status due to surgical cancellations.The rest of the 18 patients were categorized into three pathological response groups, including those with pathological complete response (pCR), major pathological response (MPR), and incomplete pathological response (IPR) (Fig. 1A; Additional file 1: Fig. S1A).The clinical characteristics of the remaining 18 patients with a definitive response to NAT are shown in Additional file 3: Table S1.

Single-cell gene expression atlas of ESCC
To describe the transcriptional expression atlas of ESCC at the single-cell level, we performed scRNA-seq on all 46 samples, including tumor and normal adjacent tissues before and after NAT (Fig. 1A).By comparing the expression of canonical cell-type marker genes, we identified 14 distinct cell subsets totaling over 250,000 cells, including cancer, immune (including T, B, plasma, mononuclear phagocytes (MP), erythrocyte and mast cells), stromal (including endothelial cells (EC), mural cells and fibroblasts), parenchymal (including basal, gland mucous cells, and gland ductal cells), and Schwann cells (Fig. 1B-D; Additional file 1: Fig. S1).In post-treatment ESCC tumors, we observed a notable reduction in the proportion of cancer cells, accompanied by an augmentation in the presence of T cells (Fig. 1B).These results indicate a dynamic change within the ESCC TME in response to treatment.Meanwhile, UMAP embeddings of pre-treatment ESCC tumors, stratified by pathological response, displayed extensive overlap, implicating a profoundly intricate and heterogeneous TME among these patients (Fig. 1E).

Plasticity of ESCC epithelial cells and their developmental trajectories into cancer cells
To analyze the developmental process and the expression heterogeneity of the malignant compartment, cancer cells were isolated based on the expression of canonical marker genes from ESCC tumors (Fig. 2A; Additional file 1: Fig. S2A).Through unsupervised clustering analysis, twelve distinct clusters of cancer cells were identified, and the extent of chromosomal copy number variations (CNVs) across the entire genome was inferred through a comprehensive analysis of scRNA-seq data and calculation of CNV scores [54,55] (Additional file 1: Fig. S2B,  C).Compared to non-malignant cells, all cancer cell clusters displayed heterogeneity in gene expression and CNV status.Notably, clusters 5 and 10 showed higher CNV scores compared to others, indicating a greater degree of chromosomal instability in their genomic profiles.The dynamic characteristics and heterogeneity of malignant cancer cells were further investigated by trajectory analysis.To determine the start point of the inferred pseudotemporal cancer cell ordering, we first applied SLICE (Single Cell Lineage Inference Using Cell Expression Similarity and Entropy) to quantify singlecell entropy (scEntropy) of each cancer cell type [44].As entropy inversely correlates with cell differentiation state, cells with lower scEntropy had more well-defined cell fates and functionalities compared to those clusters with higher scEntropy.In this context, clusters 1 and 11 exhibited a significantly lower entropy compared to others, marking the start point in trajectory mapping (Additional file 1: Fig. S2D).Subsequently, we used the "Monocle 2" method [56][57][58][59][60] to construct entropy-directed transitional trajectories of cancer cells (Fig. 2B).Notably, differentiation trajectory analysis revealed a branched structure of cancer cells, with several clusters of cancer cells positioned along the extended branch away from the start point, suggesting that these cell clusters were likely associated with higher functional uncertainty and greater differentiation potential (Fig. 2C).Interestingly, we observed that less differentiated clusters (higher entropy, such as clusters 2, 4, 5, 9, 10) were primarily identified in pre-treatment tumor samples of pCR and MPR patients, while more differentiated clusters (lower entropy, such as clusters 1 and 11) were predominantly found in posttreatment tumor samples (Fig. 2D).These observations suggest an interesting link between cancer cell differentiation and the pathological response to neoadjuvant chemo-immunotherapy.Specifically, ESCC tumors enriched with less differentiated cancer cells might be associated with a better pathological response to neoadjuvant chemo-immunotherapy, whereas tumors enriched with clusters of more differentiated cancer cells may resist the treatment.Likewise, the pseudotime faceted map of paired tumors categorized by response displayed a consistent pattern, with a better pathological response observed in tumors enriched with cancer clusters characterized by higher entropies (Additional file 1: Fig. S2E).Nevertheless, it is crucial to recognize that scRNA-seq inherently presents limitations in precisely calculating cell proportions, as demonstrated by the underrepresentation of cancer cells in post-treatment tumor samples from IPR patients.Consequently, these results should be approached with caution, and future validation is imperative.
Furthermore, we analyzed the transcriptomic changes of cancer cells associated with transitional states (Fig. 2E,  F).Three phases of differentiation were identified.Cancer cells in state 3, which positively responded to NAT, exhibited significantly elevated expression levels of immune-related genes, including those implicated in neutrophil activation, antigen processing and presentation, and response to IFNγ (Fig. 2F).Conversely, cancer cells mapped to state 2 showed significantly increased expression of genes associated with mRNA catabolic processes.On the other hand, state 1 was predominantly occupied by cancer cells characterized by higher expression of genes involved in epidermal development.These results were consistent with the dynamic changes in gene expression of cancer cells during the transition (Additional file 1: Fig. S2F, G).Collectively, these findings suggest that state 3 likely encompassed cancer cells actively engaged in immune responses, whereas state 1 appeared to consist primarily of epithelial-like cells at the transient stage.

Dissecting the transcriptomic inter-tumor heterogeneity of cancer cells in ESCC
Hotspot analysis [46] was performed to explore informative gene modules characterized by distinct expression patterns and their correlation with pathological response to NAT.It is noteworthy that unsupervised clustering of cancer cells, as a data-driven method, facilitates the identification of subpopulations of cancer cells sharing similar expression patterns.On the other hand, hotspot analysis represents a hypothesis-driven approach, wherein specific gene sets and pathways are predefined based on prior knowledge or biological significance.Here, we identified 12 distinct gene modules corresponding to the 12 previously annotated cancer cell clusters in ESCC (Fig. 2G; Additional file 1: Fig. S2H).Cancer cluster 3 was predominantly observed in pre-treatment tumor samples of IPR patients and exhibited a strong correlation with module 3 (Fig. 2D; Additional file 1: Fig. S2H).However, it remains inconclusive whether the significant reduction of cluster 3 in post-treatment samples is correlated with the efficacy of NAT due to the limited representation of cancer cells in IPR_T_A samples (Fig. 1E).On the other hand, cancer cluster 1 with lower entropies was significantly correlated with module 1 (Jaccard similarity coefficient = 0.4).Consistent with the notion that cluster 1 cells displayed a stronger resemblance to normal epithelial cells, Gene Ontology (GO) enrichment analysis revealed that module 1 was dominated by genes associated with epidermal cell differentiation and development (Additional file 1: Fig. S2I).Notable genes within this category included EMP1, SCEL, SPINK5, ANXA1, SPRR3, PPL, CNFN, and TGM3.In contrast, gene modules associated with pathological complete response to NAT were primarily implicated in biological processes such as extracellular matrix organization (module 6), and activation of protein complex assembly (module 9) (Additional file 1: Fig. S2J).Regarding gene modules significantly contributing to MPR, the results of GO functional analysis unveiled a noteworthy enrichment of genes associated with epidermis development and response to type I interferon (module 7), response to hypoxia (module 10), and protein targeting to ER (module 4) (Additional file 1: Fig. S2K).
Given the unfortunate unavailability of survival data from the study cohort, we conducted an exploratory analysis utilizing the TCGA-ESCA dataset (N = 159).A Cox regression model was applied to investigate the relationship between cancer-associated gene modules and the patient's overall survival (OS).Remarkably, higher expression of gene modules 6 and 12 was significantly associated with a more favorable prognosis of patients (P = 0.039 and 0.044, respectively; Additional file 1: Fig. S2L, M).In contrast, cluster 10-associated module 9 and cluster 5-associated module 10 were likely associated with a worse prognosis, though the difference did not reach statistical significance (Additional file 1: Fig. S2H,  L).Nevertheless, it is worth noting that bulk RNA-seq data was used rather than scRNA-seq data and that the external dataset adheres to a different treatment regimen from the clinical settings in our analysis.Further exploration and validation of the association between cancer cell differentiation and immunotherapy-related patient prognosis are warranted in future studies.
We also employed SCENIC [49] to identify specific combinations of transcription factors that govern the expression of their respective target genes within cancer cells at the single-cell level.The top 5 regulons in cancer cells from patients exhibiting various responses to neoadjuvant anti-PD-1 combination therapy showed substantial differences in their expression patterns (Fig. 2H).Notably, pre-treatment pCR tumors (pCR_T_B) showed an upward trend in the expression of ZNF217 and MEF2A, whereas a downward trend was observed in the expression of TFCP2L1, ZSCAN31, and SOX2 compared to the other two groups (all P < 0.001, Jonckheere-Terpstra test) (Fig. 2I).Importantly, SOX2 has been extensively studied in the context of cancer, playing a pivotal role in tumor cell proliferation, cancer metastasis, and response to therapies [61].Accumulating evidence suggests that SOX2 promotes cancer cell stemness, and its overexpression is associated with lymph node metastasis and poor prognosis in cancer patients [62].Collectively, these findings suggest that the dynamic change in these transcription factors and their associated regulatory networks might be associated with the pathological response of patients who received NAT.
While PD-L1 expression has been employed to predict the efficacy of pembrolizumab in advanced EC [63], the predictive role of PD-L1 expression remains inconclusive for neoadjuvant therapy in EC patients.In our study, a total of 14 patients exhibited TPS scores < 50% or CPS scores < 50, thereby being classified as having PD-L1 low expression.Among these patients, 9 individuals (64.3%) exhibited pCR or MPR to NAT and five demonstrated IPR examined by conventional histologic approach (Additional file 3: Table S1).Dissecting the transcriptomes of cancer cells at single-cell resolution in low PD-L1 patients may provide significant clinical insights for understanding how these patients respond to neoadjuvant chemo-immunotherapy.In patients exhibiting low PD-L1 expression, a notable upregulation of genes associated with the protein targeting pathway was observed in tumors that achieved pCR.Conversely, in IPR tumors, there was an increased expression of genes associated with responding to oxidative stress (Additional file 1: Fig. S2N-P).These findings collectively suggest a plausible connection between the activation of these transcriptional signatures in cancer cells and the different pathological responses to NAT observed in patients with low PD-L1 expression.

Increased immune surveillance associated with better pathological response revealed by ESCC TME
We analyzed cell-cell communication across all cell types within the ESCC TME using the "CellphoneDB" method.The results revealed a positive correlation between TME activity and the patient's pathological response to NAT, with a notably higher number of ligand-receptor pairs (LRPs) observed in baseline pCR tumors (Fig. 3A).We then categorized T cells based on canonical marker gene expression, identifying nine subsets totaling 7209 cells (Fig. 3B, C).Particularly noteworthy was the significantly greater proportion of CD8 + effector T cells in baseline pCR tumors compared to IPR tumors (Fig. 3B).Meanwhile, T cells exhibited stronger interactions with both MP cells and within their own population (Additional file 1: Fig. S3A, B).These findings suggest a positive association between CD8 + effector T cell infiltration and improved pathological response to NAT in ESCC.

Elevated expression of IFNγ-responsive gene signature in CD8 + effector T cells predicts pathological response of ESCC patients to NAT
Transcriptional analysis was performed to identify treatment response-related molecular signatures in CD8 + effector T cells at the single-cell level.Emphasis was placed on identifying DEGs and enriched pathways within pre-treatment ESCC tumors of pCR and IPR patients who exhibit the most pronounced variability in response to NAT (Fig. 3D).The expression levels of genes involved in the cellular response to IFNγ (such as MT2A, CCL4, and B2M) were notably upregulated in baseline pCR tumors compared to those in IPR tumors (Fig. 3E, F).Of note, a decreased expression of metallothionein (MT) family genes, including MT2A, MT1X, and MT1E, has been previously associated with resistance to anti-PD-1 therapy in ESCC [64].While MT2A expression was significantly higher in MPR patients compared to those with IPR (P = 0.008), the difference did not reach statistical significance (Additional file 1: Fig. S3C-F).These findings suggest that the functional significance of CD8 + effector T cells in baseline MPR patients may resemble that in pCR patients rather than IPR patients.Assessing DEGs across various comparison groups, we found that upregulated genes linked to better pathological responses were mainly markers associated with cellular response to IFNγ (CCL4, MT2A, and B2M), neutrophil degranulation (B2M and HSPA8), regulation of I-kappaB kinase/NF-kappaB signaling (S100A4), and negative regulation of T cell apoptotic process (TSC22D3) (Fig. 3G).In contrast, downregulated gene signatures in baseline ESCC tumors of patients with better NAT-associated pathological response were primarily involved in the immune response-activating cell surface receptor signaling pathway (such as IGHG4, IGHA1, C3AR1, FCER1G, TYROBP, AAPL1, and LPXN) (Fig. 3E; Additional file 1: Fig. S3G).Among these genes, higher expression of C3AR1, FCER1G, and AAPL1 has been previously linked to an unfavorable prognosis in ESCC patients [65,66].
To gain further insight into potential treatment strategies, we examined the top 30 ligand-receptor interactions between CD8 + effector T cells and cancer cells, as well as other T cell subsets.We observed a significantly lower number of cell-cell interactions in baseline tumors of IPR patients compared to those with pCR (Fig. 3H,  I).In both IPR and pCR tumors, CCL5 was commonly expressed as a ligand on CD8 + effector T cells, with its receptor CCR1 expressed on NK cells and Treg cells.However, the identification of CCL5-CCR1 interactions within CD8 + effector T cells was observed only in baseline ESCC tumors of pCR patients.In contrast, CCL5-CCR5 interactions were widely expressed across various cell types in pCR tumors at baseline, except for cancer cells, NK cells, and innate lymphoid cells (ILC).Besides, CCL3 and CCL4 may also stimulate cell-cell communications through most immune cell subsets by binding to its receptor CCR5.Through assessment of the diverse expression patterns of immune checkpoints, we noticed a significant enrichment of interactions, including TIGIT-PVR and LGALS9-HAVCR2 in baseline IPR tumors, suggesting novel immunotherapy targets for ESCC patients who may not benefit from NAT (Fig. 3J, K).
Next, transcriptional changes in CD8 + effector T cells associated with pathological response to NAT were examined using pCR and IPR patients with low PD-L1 expression.Consistently, pCR patients with low PD-L1 expression exhibited higher expression levels of genes involved in the cellular response to IFNγ, along with lower expression of genes associated with immune response-activating signaling transduction pathway (Additional file 1: Fig. S3H-J).Overall, the recapitulation of GO enrichment analysis in low PD-L1 patients suggested that increased expression of IFNγ-responsive gene signatures and decreased expression of genes related to immune response-activating signaling pathway in CD8 + effector T cells might predict a favorable pathological response in locally advanced ESCC patients who underwent NAT.

NK cells constitute a minor population but induce a cytotoxic cellular response in the ESCC TME
NK cells constitute a unique subset of lymphocytes larger than T cells and B cells, capable of rapidly targeting and eliminating specific cells without prior immunization or MHC restriction [67,68].Our study identified a lower proportion of NK cells in pCR patients compared to those in non-pCR patients following NAT (Fig. 3B).Similarly, we assessed DEGs in the subgroup analysis (Fig. 4A,  B; Additional file 1: Fig. S4A-C).Pathway enrichment analysis revealed that genes significantly upregulated in patients exhibiting a positive pathological response to NAT were primarily involved in enhancing lymphocyte chemotaxis (such as CCL4, CCL3, and XCL2) and regulating T cell-mediated immunity (B2M and KLRD1) (Additional file 1: Fig. S4A).Compared to IPR patients, those achieving pCR showed downregulation in genes associated with neutrophil-mediated immunity pathways, including CXCL1 and ITGB2 (Fig. 4B).
We then investigated the interactions between NK cells and various cell types.The expression of LRPs involving NK cells and immune cells, such as CCL5-CCR4, CCL5-CCR1, and CCL4-CCR8, were exclusively observed in pCR patients before NAT but not in IPR patients (Fig. 4C, D).Furthermore, while the CCL5-CCR1 interaction was detected between CD8 + effector T cells and NK cells in both pre-treatment pCR patients and IPR patients, this interaction pair was also identified in pCR patients between CD8 + naïve T cells and NK cells, as well as within the NK cell population.A range of immune checkpoint pairs, such as TNFSF4-TNFRSF4, CTLA4-CD86, LGALS9-HAVCR2, and TIGIT-PVR, between NK cells and CD8 + effector/exhausted T cells, were significantly enriched in pre-treatment IPR patients, indicating additional potential targets for immunotherapy in nonresponders following the current treatment regimen (Fig. 4E, F).
Similar to what was found in CD8 + effector T cells, in low PD-L1 expression pCR patients, NK cells exhibited a markedly elevated expression of genes involved in the cellular response to IFNγ, as well as those implicated in the response to interleukin-1, such as CCL4, CCL3, CCL5, XCL1, UBB, ANXA1, and XCL2 (Additional file 1: Fig. S4D).Conversely, in IPR patients with low PD-L1 expression, there was a significant upregulation of genes associated with antiviral responses and those pivotal in RNA catabolic processes.

Decreased proportions of Treg cells in ESCC indicate better pathological response following NAT
Treg cells, originating in the thymus, exert regulatory controls over various immune cells including T cells, B cells, NK cells, dendritic cells (DCs), and macrophages through both humoral and cell-cell contact mechanisms [69].The infiltration of Treg cells is linked to poor prognosis due to their fundamental suppression mechanisms involving CTLA4, IL-2, IL-10, and other immune-suppressive cytokines and substances [69,70].In this study, we observed a decline in Treg cell numbers among ESCC patients exhibiting better pathological responses to NAT, consistent with their suppressive role within the TME (Fig. 3B).Then, we performed DEG analysis between patient subgroups categorized based on their pathological response to NAT (Fig. 5A, B; Additional file 1: Fig. S5).Comparison of baseline tumor samples from the two patient cohorts with the most variable responses to NAT revealed a significant downregulation of genes in pCR patients involved the humoral immune response (such as S100A8, IGHG4, S100A9, IGKC, IGLC3, IGHG3, RGCC , HLA-A, IGHG1, CXCL8, CXCL1, A2M, KRT6A, IGHA1, and CXCL2), B cell activation (such as IGHG4, IGKC, IGLC3, IGHG3, IGHG1, BATF, TNFRSF4, RASGRP1, IGHA1, NDFIP1, ZFP36L1, TBC1D10C, andCTLA4), and myeloid cell differentiation (such as JUN, HSPA1B, HSPA1A, JUNB, ZFP36, ADAR, PIP4K2A, NFKBIA, NR4A3, FOS, and GNAS) (Fig. 5B).As protein targeting, ribonucleoprotein complex biogenesis, and mRNA catabolic process were identified as the top three signaling pathways significantly enriched in pCR patients, we reasoned that Treg cells may contribute to improved clinical responses to NAT by maintaining cellular homeostasis, fostering self-tolerance, and facilitating tissue repair.
The examination of cell-cell interactions revealed that Treg cells displayed restricted interactions with other T cell subsets in IPR patients, whereas the communication was slightly intensified in pCR patients by the presence of CCL5-CCR1 and CCL4-CCR8 interactions (Fig. 5C, D).Immune checkpoint pairs, such as TIGIT-PVR, TNFSF4-TNFRSF4, LGALS9-HAVCR2, and CD70-CD27, were identified in baseline IPR patients as potential targets for immunotherapy (Fig. 5E, F).

Discussion
Recent studies have undertaken scRNA-seq analyses to uncover the mechanisms driving the response of esophageal cancer to neoadjuvant chemotherapy with or without concurrent radiotherapy contributed by immune cells in the TME [30,31].Nevertheless, it is crucial to emphasize that nearly half of patients still encounter local recurrence or distant metastases after preoperative neoadjuvant chemoradiotherapy followed by surgery [20].Hence, exploring novel and effective treatments remains an urgent need for improved survival of EC patients.Our study employed scRNA-seq to evaluate immune contexture and their correlation with pathological response to the emerging neoadjuvant chemo-immunotherapy, which holds significant clinical implications for the precise treatment of resectable ESCC patients.
The present investigation has yielded several significant findings.Firstly, we have deciphered the complicated compositions of the TME in ESCC using scRNA-seq, which identified 14 major cell clusters, including cancer, immune, stromal, and Schwann cells.At baseline, the TME was predominantly composed of epithelial cells, particularly cancer cells, followed by MP cells.Neoadjuvant chemo-immunotherapy may modulate cell type contexture in several ways.Most notable was a reduction in the proportion of cancer cells and an increase in T cells following therapy, indicating a dynamic change within the ESCC TME in response to treatment.Previous research has shown that neoadjuvant chemotherapy, whether administered alone or in combination with neoadjuvant radiotherapy, resulted in an augmentation of CD8 + T effector cell infiltration and a reduction in Treg cell populations in esophageal cancer patients [30,31].Consistently, we observed similar modulations of these T cell repertoire stratified before or after NAT, suggesting a potential synergy of different types of preoperative neoadjuvant therapies.Furthermore, patients across different response groups showed both intra-and inter-tumoral heterogeneity in their pre-treatment tumor samples, emphasizing the need for a more comprehensive investigation of the heterogeneous responses of malignant and immune cell components to NAT.Nonetheless, it is important to acknowledge that scRNA-seq has inherent limitations in accurately determining cell proportions.Therefore, conclusions drawn regarding proportions should be approached with caution.As evidence, we observed an underestimated proportion of cancer cells in post-treatment tumor samples from IPR patients.
Secondly, we discovered a potential link between tumor cell differentiation states and various major expression programs, which may potentially impact the pathological response to neoadjuvant chemo-immunotherapy.Worth noting that single-cell trajectories were constructed using the "Monocle 2" algorithm to elucidate the diverse cell fates of cancer cells within ESCC tumors.Our study did not propose any direct evolutionary/developmental relationships among individual tumors that are inherently genetically distinct.The results showed that the more differentiated cancer cells located at the pre-branch (state 1) may exhibit a closer resemblance to normal epithelial cells.On the contrary, less differentiated clusters (higher tumor stemness) positioned at the terminal point in the trajectory may display a stronger likeness to stem cells, characterized by a higher level of functional uncertainty.These cells were likely associated with neutrophil-mediated immunity and the processing and presentation of antigens.Notably, they are more prone to exhibiting a favorable pathological response to NAT.In particular, we observed a remarkable decrease in the number of cells within post-NAT tumors in cancer clusters 2, 4, 5, 9, and 10, whereas clusters 1 and 11 may potentially represent cell subsets demonstrating an adaptive immune resistance to NAT.
Thirdly, we utilized specific immune cell populations and performed a comprehensive analysis of DEGs and pathways among subgroup patients manifesting different pathological responses to NAT, with particular emphasis on pCR and IPR patients who demonstrate the most divergent clinical responses to therapy.Notably, Carroll et al. have reported an ICI-responsive gene signature known as INCITE [24].The upregulation of INCITE genes involved in various inflammatory and IFNγ-related pathways could induce tumor shrinkage in EAC patients.Consistent with this finding, we noted a significant upregulation of genes associated with cellular response to IFNγ in both CD8 + effector T cells and NK cells of baseline pCR tumors compared to that in IPR tumors.Among these DEGs, MT2A may serve as a promising favorable prognostic indicator in ESCC patients undergoing neoadjuvant anti-PD-1 combined therapy, whereas C3AR1 [65], FCER1G, APPL1 [66], S100A8 [71], S100A9, CXCL1 [72], and ITGB2 [73] might be associated with an unfavorable prognosis in ESCC.The finding that IFNγ-responsive genes correlate with better pathological response to NAT suggests that the immune response plays a crucial role in determining treatment efficacy.The potential of IFNγ to stimulate lymphocyte and macrophage functions highlights its significance as a therapeutic target in treating ESCC patients.
Furthermore, T cells express a broad spectrum of ligands and receptors for chemokines, cytokines, checkpoint molecules, and growth factors.The aberrant expression of the C-C chemokine ligand 5/C-C chemokine receptor type 5 (CCL5-CCR5) axis has been recognized as a pivotal contributor to tumor progression, facilitating a more conducive TME for hematological malignancies and various solid tumors [74,75].However, CCL5 may also augment the immunotherapy response by recruiting antitumor T cells and dendritic cells to the TME [76][77][78][79].Here, we demonstrated extensive interaction between CCL5 and its specific receptor CCR5 in various immune cells of pre-treatment pCR patients.In contrast, very low levels or absence of the CCL5-CCR5 interaction pair could be detected in baseline tumor samples of IPR patients.These results suggest that the binding of CCL5 to its specific receptor CCR5 may facilitate cooperation among diverse immune cell types, stimulating the cytotoxic function of CD8 + effector T cells.It is also worth mentioning that CCR5 is a promiscuous G protein-coupled receptor that binds with high affinity to not only CCL5 but also CCL3 and CCL4 [80].Consistently, CCL3-CCR5 and CCL4-CCR5 interaction pairs were also detected between CD8 + effector T cells and various T cell subsets in our study.Interestingly, we noted a correlation between higher CCL4 expression in CD8 + effector T cells and a more favorable pathological response to NAT, even in the comparison between pCR and MPR groups (Fig. 3G).This finding aligns with a prior study demonstrating a link between elevated CCL4 expression, increased intra-tumoral CD8 + T cells, and prolonged overall survival in EC tumors, likely due to the capacity of CCL4 to attract CCR5 + cytolytic lymphocytes [81].Besides, the presence of CCL4 in NK cells, along with its receptor SLC7A1 in cancer cells, and CCR8 in CD4 + naïve T cells/Treg cells signified an activation of NK cells [82], which may predict a more favorable pathological response in pCR patients before NAT.
Despite neoadjuvant anti-PD-1 combination therapy showing promising clinical benefits in reducing tumor cell contents, a substantial subset of patients exhibited minimal response to the treatment and were categorized as IPR patients in our study.We reported a series of immune checkpoint interaction pairs, including CTLA4-CD86, TIGIT-PVR, LGALS9-HAVCR2, and TNFSF4-TNFRSF4, which could offer further insights for the discovery and development of new therapeutic agents for ESCC.
While our study provided an unbiased and highprecision scRNA-seq analysis to decipher the TME in ESCC patients undergoing neoadjuvant chemo-immunotherapy, it is crucial to acknowledge the limitations of this pioneering study.One notable limitation was the restricted sample size, especially the scarcity of adequate paired samples for each patient before and after NAT, potentially limiting our ability to discern differences among various pathological response subgroups in ESCC patients.However, it is important to recognize the challenges associated with collecting paired samples for single-cell analysis, considering the relatively innovative nature of the treatment approach.Clinical constraints, such as the necessity for patients to receive the same treatment, the requirement for samples to meet sequencing criteria, and the technical demands of scRNA-seq necessitating fresh sample dissociation, all add complexity to the process.Nevertheless, it is imperative to address this limitation by expanding the sample size in future investigations.Furthermore, the absence of external validation through techniques such as immunohistochemistry or immunofluorescence, or comparative analysis with an external dataset involving patients undergoing the same treatment, might potentially impact the generalizability and overall validity of our findings.Therefore, our findings remain exploratory and should be approached with caution, considering the limited sample size and the extent of validation provided.Further exploration and validation of our results are warranted in future studies.

Conclusions
In conclusion, our study provides a comprehensive single-cell atlas of ESCC based on large-scale scRNA-seq results.We thoroughly assessed the development of cancer cells and T cell transcriptome, revealing promising biomarkers associated with the pathological response to NAT.Our findings have significant implications for drug development and precision medication in ESCC patients.

Fig. 1
Fig.1 Single-cell atlas of ESCC.A Schematic demonstration of the experimental workflow for single-cell RNA sequencing and computational analysis.B UMAP embedding overlaid with unsupervised cluster cell type annotations (left), and sample origin annotations (right).Pie charts (top) demonstrate the proportion of each cell cluster based on the cell type or sample origin classifications.C Average expression profile of canonical marker genes to separate cell populations across 46 samples collected from 22 ESCC patients.D The proportion of cell subsets across 46 samples (left) and the proportion change of cells across 12 sample groups (right).E UMAP profiles delineating cell subsets collected from ESCC tumors, categorized based on treatment status and pathological response to neoadjuvant chemo-immunotherapy. pCR, pathological complete response; MPR, major pathological response; IPR, incomplete pathological response; T, tumor; N, normal tissue; B, pre-treatment; A, post-treatment

Fig. 2
Fig. 2 Transitional states of cancer cells revealed by trajectory mapping.A UMAP embedding of cancer cells overlaid with unsupervised cluster cell type annotations (left), proportional sample contributions to each cell type cluster (middle), and sample label (right).B Pseudotime trajectory of ESCC cancer cells in a two-dimensional state space inferred by the "Monocle 2" method.C Cancer cells mapped to the branched structure in the trajectories (left) and the distribution of cancer clusters in ESCC tumors stratified by treatment and pathological response (right).D Cell number count in cancer clusters identified in ESCC tumors of patients before (T_B) and after (T_A) neoadjuvant chemo-immunotherapy.E Pseudotime trajectory of ESCC cancer cells showing three differentiation states.F Investigation of biological processes (BP) through Gene Ontology (GO) pathway enrichment analysis in cancer cells across three differentiation states.G Heatmap shows the top 500 genes with significant autocorrelation grouped into 12 gene modules based on pairwise correlations of gene expression in cancer cells.H The scatter plots show the expression of the top five highly expressed regulons in each of the six cell subsets.I Heatmap of the area under the curve (AUC) scores of expression regulation by transcription factors in ESCC tumors stratified by treatment and pathological response, estimated by SCENIC

Fig. 4 Fig. 5
Fig. 4 Transcriptional analysis of NK cells in baseline tumors of pCR and IPR patients.A Differentially expressed genes (DEGs) in subgroup analysis.Significant P values were labeled in red, and the y-axis represents the log2 fold change.B GO and KEGG pathway enrichment analysis for baseline ESCC tumor in pCR patients compared to IPR patients.The length and color of the bars represent enrichment significance and classifications.The number inside the spheres represents the number of related mRNAs enriched in the specific pathway.C, D Dot plots show the top 30 ligand-receptor interactions from NK cells obtained from baseline ESCC tumors of IPR (C) and pCR patients (D).The size of the circle represents the P values, and the color of the circle indicates the average expression level of interacting pairs.The cell clusters labeled in blue and red on the x-axis indicate that NK cells act as ligands and receptors in the interaction pairs, respectively.E, F Dot plots show the immune checkpoint ligand-receptor interactions from NK cells obtained from baseline ESCC tumors of IPR (E) and pCR patients (F)

Table 1
Clinical characteristics of patients (N = 22) IPR incomplete pathological response, MPR major pathological response, pCR pathological complete response, ESCC esophageal squamous cell carcinoma, EASC esophageal adenosquamous carcinoma, TPS tumor proportion score, CPS combined positive score